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Abstract 

Background: The diffusion and reaction of the transmitter acetylcholine in neuromuscular junctions and the 
diffusion and binding of Ca 2+ in the dyadic clefts of ventricular myocytes have been extensively modeled by 
Monte Carlo simulations and by finite-difference and finite-element solutions. However, an analytical solution that 
can serve as a benchmark for testing these numerical methods has been lacking. 

Result: Here we present an analytical solution to a model for the diffusion and reaction of acetylcholine in a 
neuromuscular junction and for the diffusion and binding of Ca 2+ in a dyadic cleft. Our model is similar to those 
previously solved numerically and our results are also qualitatively similar. 

Conclusion: The analytical solution provides a unique benchmark for testing numerical methods and potentially 
provides a new avenue for modeling biochemical transport. 



1. Background 

In intercellular and intracellular spaces, passive transport 
of biomolecules is a common phenomenon. Because 
such processes are difficult to probe directly by experi- 
ments, numerical modeling is increasingly used to gain 
insight. Two processes that have been extensively mod- 
eled are the diffusion and reaction of the transmitter 
acetylcholine in a neuromuscular junction [1-6] and the 
diffusion and binding of Ca + in the dyadic cleft of a ven- 
tricular myocyte [7,8]. In contrast to previous numerical 
approaches, here we present an analytical solution of a 
model for the diffusion and reaction of acetylcholine in a 
synaptic cleft (or Ca + in a dyadic cleft). Our model is 
similar to those previously solved numerically; hence our 
analytical solution potentially provides a new avenue for 
modeling biochemical transport. More importantly, an 
analytical solution provides a unique benchmark for test- 
ing numerical methods. Such a solution has been lacking 
up to now; the present work fills this gap. 

Neuromuscular junction refers to the cleft between a 
motor neuron and a muscle fiber. As illustrated in 
Figure 1, the neuronal signal for muscle contraction is 
mediated by acetylcholine. These neurotransmitter 
molecules are initially inside vesicles located in the pre- 
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synaptic axon terminal. When an action potential 
reaches the axon terminal, the vesicles release acetylcho- 
line molecules into the synaptic cleft. These molecules 
then diffuse toward the post-synaptic membrane and 
bind to acetylcholine receptors in the membrane. Acet- 
ylcholine binding activates these ligand-gated ion chan- 
nels, allowing Na + to flow in and generating an action 
potential along the muscle fiber. Finally excess acetyl- 
choline molecules around the post-synaptic membrane 
are broken down by acetylcholinesterase to prevent con- 
tinued activation of acetylcholine receptors. 

A related system is a dyadic cleft, which spans the gap 
between the cell membrane in a transverse tubule and 
the membrane of a sarcoplasmic reticulum. As Figure 2 
illustrates, Ca 2+ can enter the cell through L-type Ca 2+ 
channels on the cell membrane in response to the arri- 
val of an action potential. The ions then diffuse to reach 
and activate ryaonodine receptors in the membrane of 
the sarcoplasmic reticulum. The activated ryaonodine 
receptors release Ca + from the sarcoplasmic reticulum, 
which ultimately lead to muscle contraction. 

Here we propose a simple but not unrealistic model for 
the diffusion and reaction of acetylcholine in a neuromus- 
cular junction. The model also applies to the diffusion and 
binding of Ca + in a dyadic cleft. We are able to find an 
analytical solution for this model. For convenience we 
describe our model in the language of neuromuscular 
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Figure 1 Illustration of a neuromuscular junction. Presented are the key players in the diffusion and reaction of the neurotransmitter, 
acetylcholine (ACh). 



junction. As shown in Figure 3, we model the synaptic 
cleft as the space between two infinite, parallel planar 
membranes. In the pre-synaptic membrane, there is a peri- 
odic array of circular openings, via which acetylcholine 
molecules enter the cleft. In the post-synaptic membrane, 
there is a periodic array of circular disks, where the acetyl- 
choline molecules are absorbed. The quantity of interest is 
the total flux, J(t), at time t of acetylcholine molecules 
across the post-synaptic membrane. This model allows us 
to use periodic boundary conditions in the transverse 
direction. Our results are qualitatively similar to those 



obtained previously by numerical solutions [1-7]. More 
realistic ingredients can be added to our model and still 
permit analytical solutions. 

2. Methods 

We set up a coordinate system such that the x and y 
axes are parallel to the pre- and post-synaptic mem- 
branes. The synaptic junction has depth L z . The junc- 
tion is periodic in the x and y directions, with 
periodicities of L x and L y , respectively. In each "unit 
cell", a synaptic vesicle bursts at time t = 0, releasing 
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Figure 2 Illustration of a dyadic cleft. Presented are the key players in the diffusion and binding of Ca 
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Figure 3 Our model for both a neuromuscular junction and a dyadic cleft. Panel A: The pre-cleft membrane (top) contains a periodic array 
of disks for the influx of ligands; the post-cleft membrane (bottom) contains a periodic array of absorbing disks representing receptors. Panel B: 
The dimensions of a unit cell. 



the neurotransmitters into the cleft. We model the 
release of the neurotransmitters as a transient flux, u(t), 
that is confined to a circular opening with radius R. We 
place the synaptic vesicle at the center of the presynap- 
tic face of the unit cell. After diffusing to the post- 
synaptic membrane, neurotransmitters are absorbed by 
a circular disk in each unit cell, with radius a. We place 
this "sink" also at the center of the post-synaptic face of 
the unit cell. The exact shapes and locations of the pre- 
synaptic opening and the post-synaptic sink are not 
essential for the analytical solution of our model. The 
quantity of interest is the total flux, J{t), through the 
post-synaptic face of each unit cell. 

We place the origin of the Cartesian coordinate sys- 
tem at the center of the pre-synaptic face, with the z 
axis pointing toward the post-synaptic face. The concen- 
tration of neurotransmitters at position r and at time t 
is C(r, t). Within the synaptic junction, it is governed by 
the diffusion equation, 

e> C M) = DV2 C (r,0 (1) 
dt 

where D is the diffusion constant. The initial condi- 
tion is 

C(r, 0) = 0 (2) 
The boundary condition at z = 0 is 

D ^ = u[t) when p < R (3a) 
dz 

D ac(r,t) =Q whenp> £ (3b) 

dz 

where p = (x 2 + y 2 ) 112 is the distance to the z axis. 
The boundary condition at z = L z is 



C(r, t) = 0 when p < a (4a) 

p dC(r,t) =Q whenp>a (4b ) 
dz 

We solve the problem in Lapace space. For a function 
f(t) of t, we denote the Laplace transform as 

f{s) = jdte~ st f(t) . The Laplace transform of Eq. (1), 

o 

using the initial condition of Eq. (2), is 

sC{t,s) = DV 2 C{t,s) (5) 

The solution appropriate for the periodic boundary 
conditions in the x and y directions has the form 

C(r, 5 ) = £ [a ta ( S >'-W« + PU^ r " Mz ]co S (2fcx/ L x )cos(2mny / L y ) (6) 
i,m=0 

where 

y lm (s) = [ s I D + {2h I L,) 2 + (imny / L y ) 2 (7) 

The coefficients Ui m {s) and Pi m (s) are to be determined 
from boundary conditions. 

To make use of the boundary condition of Eq. (3), we 
need the 2-dimensional cosine transform of a function in 
x and y that has value 1 when p < R and value 0 when 
p > R. The coefficient of the cos(2ln/L x )cos(2mn/L y ) term 
is 

<?/ m = jj dxdy cos{2lnx / L x ) cos(2mny / L y ) ^ 

X y p<R 

where e 0 = 1/2 and £/ = 1 for / > 0. For later reference, 
we note that the function having value 1 when p < a 
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and value 0 when p > a has the following coefficients in 
its 2-dimensional cosine transform: 



Pirn 



dxdy cos(2lnx / L x )cos(2mny / L y ) 



(9) 



Applying the boundary condition of Eq. (3), we find 



tin 



(10) 



For the boundary condition of Eq. (4) at z = L z , we 
use the constant- flux approximation [9]: 



D BC f r ' 5) = Q(s) when p<a 



D 



dz 

8C(r,5) 
dz 



0 when p > a 



(11a) 



(lib) 



where the quantity Q(s) is to be determined by the 
condition 



: C(r, s) {na 2 )~ l J*J dxdyC{x, s) = 0 



(12) 



p«i 



Equation (11) leads to 

DY lm (s)[a Im (sy*> ls)L * ~ p lm {s)e-^ L * ] = Q(s)p ltn (13) 

With Eqs. (10) and (13), we solve for coefficients a.i m 
(s) and Pi m {s). The results are 



Mt. _ 1 -<ji m Ks) + Pi m Q{sY y ' m(s)L * (14; 
Dy lm {s) e r„ML z _ e -y„„(s)L z 



a lm e 



a „-r.»(s)t, 1 -q lm u{s) + p lm Q{s)e n ™ {s)L > {14b) 

p,m D Ylm (s) e y^- e -y^ 

Inserting Eq. (6) in Eq. (12) and using Eqs. (14) and 
(9), we find Q(s) as 



Q(s)=u(s) 



l,m=0 



VinAlm 

(s)sinh[7, m (s)L z ] 



y p;,„ 2 cosh[7; m (5)L z ] 



(15) 



Finally the total flux through the sink on the post- 
synaptic face is 



) is) = k a 2 Q(s) 



(16) 



The flux accumulated over all times, 



is of interest. In our model, 



J = /(0) = 7iR 2 u{0) 



(17a) 



(17b) 



which is independent of a. Equation (17b) is simply a 
consequence of ligand conservation, i.e., the total num- 
ber of ligands released from the synaptic vesicle is the 
same as the total number of ligands absorbed by the 
receptors. 

The analytical solution derived above can be imple- 
mented on any function u{t) modeling neurotransmitter 
release from the synaptic vesicle. We focused on an 
exponentially decaying function: 



u(t) = e 



-t/t 0 



Its Laplace transform is 



u(s) 



1 



5 + 1/ to 



(18a) 



(18b) 



Our model then contains two parameters related to 
time: D and t 0 . For two sets of these parameters, e.g., D\ 
and t ol in one and D 2 and t 02 in the other, it can be 
shown that the corresponding response functions satisfy 



/ 1 (D 2 t) = ] 2 {D l i) when D x t ol = D 2 t c 



(19) 



We now briefly describe the details of our implemen- 
tation of the analytical solution. To calculate the qi m 
and pi m coefficients of Eqs. (8) and (9), we first carried 
out the integration over y analytically. The remaining 
integration over x was done numerically using the 
Gauss-Legendre quadrature with two points. The sum- 
mations over / and m in Eq. (15) were truncated at / = 

m = 40. The Laplace transform of /(s) was inverted by 

the Stehfest algorithm [10]. A Fortran90 code for the 
implementation is available upon request. 

3. Results 

We now present some illustrative results. The para- 
meters of our model are as follows: L x = L y = 500 nm; 
L z = 50 nm; R = 20 nm; a varied from 2.5 to 40 nm; t Q 
varied from 1 to 10 ms; and D varied in (0.4-4) x 10 
nm 2 /ms. 

Our single absorbing disk is used to model ligand 
binding to multiple receptors. Increasing a thus mimics 
an increasing number, N lec , of receptors per unit cell. 
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We expect a to increase linearly with increasing AT rec 
when N rec is small; the increase in a then slows down at 
higher A^ rec (see Discussion). Figure 4 displays the 
dependence of the response function /(f) on a (and 
hence N lec ). As a (i.e., N lec ) increases, /(f) more and 
more quickly reaches a higher and higher maximum 
and decays faster and faster. This is the expected beha- 
vior. Moreover, according to Eqs. (17b) and (18b), the 
areas under the curves for different a values are 
all equal to nR t 0 . Note also that all our /(f) curves 
have the familiar shape seen in previously numerical 
studies [1-7]. 

Figure 5 shows the change in the response function 
when the speed at which the neurotransmitters are 
released into the synaptic cleft is varied. To make a fair 
comparison, /(f) is scaled by f 0 so that, effectively, the 
total number of neurotransmitters entering the synaptic 
cleft is fixed. It is clear that, as f 0 increases (i.e., as the 
speed of neurotransmitter release decreases), the 
response function rises and decays more slowly. This 
behavior has previously been specifically modeled by 
Stiles et al. [1]. 

There is some uncertainty on the diffusion constant of 
acetylcholine in neuromuscular junctions. In previous 
models [1-6], the value of D ranged from (0.25-6) x 10 5 
nm /ms. In Figure 6 we display the change in the 
response function when D is varied from (0.4-4) x 10 
nm /ms. As expected, when neurotransmitter diffusion 
slows down, the response function is also delayed. In 
our model, decreasing D has the same effect on the 
response function as increasing f 0 [see Eq. (19)]. 

4. Discussion and Conclusion 

We have presented an analytical solution to a model 
for the diffusion and reaction of acetylcholine in a 




0 2 4 6 8 10 



t(ms) 



Figure 5 Effect of the speed of ligand release on the response 
function. Scaling of J(t) by f 0 is to ensure that the same number 
of neurotransmitters is released for all the curves with different 
f 0 values, which are shown in the figure in ms. a = 10 nm and 
D = 2 x 1 0 5 nm 2 /ms. 

neuromuscular junction. The model also applies to the 
diffusion and binding of Ca + in a dyadic cleft. Our 
results are qualitatively similar to those obtained pre- 
viously from models solved numerically [1-7]. 

Perhaps the greatest value of our analytical solution is 
that it provides a benchmark for testing numerical meth- 
ods. Diffusion and reaction of ligands in intercellular and 
intracellular spaces have been modeled either on a particle 
description or a concentration description. The former 
type of models have been solved by Monte Carlo simula- 
tions [1,7], while the latter type of models have been 
solved by either finite-difference [2,6] or finite-element 
[3-5] methods. The two types of models have been shown 
to give equivalent results [8] . The level of realism of our 
model approaches those of the models solved numerically; 
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hence our analytical solution will be able to serve as a 
good benchmark for the numerical methods. 

Our model has room for increasing the level of rea- 
lism and still allows for analytical solution. For example, 
we modeled ligand binding to receptors as absorbing. 
The binding can be modeled as partially absorbing if 
binding does not occur at every ligand-receptor encoun- 
ter. The partial absorption condition takes the form 



D d^ r ' Q = kC(i, t) when p < a 
dz 



(20) 



where the reactivity K controls the degree of partial 
absorption. When k = 0, the region p < a becomes 
reflecting and binding cannot occur. When k — > °°, the 
region p < a becomes fully absorbing and Eq. (20) 
reduces to the absorbing condition of Eq. (4a). Note 
that the bimolecular rate constant k for binding to the 
partially absorbing disk is given by [11] 



1 



1 



+ - 



1 



k ADa TtKa' 



(21) 



One can thus parameterize a and k by matching k 
with experimental data for the ligand-receptor binding 
rate constant. 

Another simplification of our model is that a single 
absorbing disk is used to represent the receptors. We 
accounted for the presence of multiple receptors per 
unit cell by increasing the radius a of this disk. One can 
identify a by requiring that the bimolecular rate con- 
stant k calculated for the single absorbing disk is the 
same as that for the N rec receptors. For the single 
absorbing disk we have k = QDa [see Eq. (21)]. If each 
of the N rec receptors is modeled as an absorbing disk 
with a small radius a 0 , then k ~ <lN rec Da 0 when AT rec is 
small [12]. Therefore a ~ N lec a 0 for small N rec . As AT rec 
increases, the rate constant k for the multiple receptors 
and correspondingly the radius a of the equivalent disk 
also increase, but the increase slows down as N rec 
becomes large [12,13]. Instead of using a single equiva- 
lent absorbing disk, analytical solution is actually still 
permitted when multiple receptors, each represented by 
a (partially) absorbing small disk, are present at arbitrary 
positions on the post-synaptic face. An alternative way 
to model the presence of multiple receptors is to assume 
that the whole post-synaptic face is partially absorbing, 
with an reactivity given by [12,13] 



K = 



N rec fe 0 

s(i-/ re J 



(22) 



where S = L x L y is the area of the post-synaptic face, 
/rec = N lec na 0 2 /S is the fraction of the post-synaptic face 
that is covered by the receptors, and k 0 is given by 



fen 



4Da 0 7nca 0 2 (l-/ rec ) 



(23) 



We have modeled ligand-receptor binding as irreversi- 
ble. This is somewhat justified for modeling the neuro- 
muscular junction, in which acetylcholinesterase can 
break down acetylcholine molecules newly released from 
the receptors. No such mechanism is present for Ca + 
in the dyadic cleft. Reversible binding can be treated by 
appropriate boundary conditions [14] on the post-synap- 
tic face. Another important detail is that both acetylcho- 
line and ryanodine receptors have multiple binding sites 
for their ligands so that there are multiple ligand-occu- 
pation states for the receptors. Again, these can be trea- 
ted by appropriate boundary conditions. 

The geometries of some of the models previously 
solved numerically are more sophisticated than that of 
our model. In particular, secondary folds of the neuro- 
muscular junction has been included in some of the 
previous models [1,3-5]. A formalism for treating 
ligand binding to a site buried in a narrow tunnel has 
been developed [15] and may be adopted for treating 
the narrow secondary folds in the neuromuscular junc- 
tion. However, analytical solution requires idealized 
geometries; the kind of realistic geometries drawn 
from electron microscopy that can be handled by 
a finite-element method [4,5] is beyond the reach of 
analytical solution. Nevertheless, with all the new 
ingredients outlined above, analytical solution will 
potentially provide a new avenue for modeling bio- 
chemical transport. 
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